/* 
 * File:   conectividad.h
 * Author: Dr. Juan Carlos Fernandez Toledano
 *
 * Created on 22 de noviembre de 2011, 18:46
 */

#ifndef CONECTIVIDAD_H
#define	CONECTIVIDAD_H

#include "objetos.h"


void conectividad(particula *part, parametros param,double rcontact, 
        int nbins, double *hist, double *cH, double *cV)
{
    /*Calcula las conectividades horizontales y verticales del sistema*/
    double dist, dx,dy,dz;
    double step = 180.0/double(nbins);
    int i,j, id;
    int count=0;
    double cosTheta, angulo, angulo2;
    /*Conectividades horizontales y verticales*/
    *cH=*cV=0.0;
    
    /*Inicializacion del histograma*/
    for(i = 0; i<nbins;i++) hist[i]=0.0;
    
    /*Recorremos todos los pares de particulas*/
    for(i=0;i<param.npart-1;i++)
    {
        for(j=i+1;j<param.npart;j++)
        {
            dist = distancia(part[i],part[j], &dx, &dy, &dz, param);
            if (dist<rcontact)
            {
                /*Las particulas estan en contacto*/
                cosTheta=dz/dist;
                /*Angulo expresado en grados*/
                angulo = 180.0*acos(cosTheta)/pi;
                angulo2 = 180.0-angulo;
                id = angulo/step;
                hist[id]+=1.0;
                id = angulo2/step;
                hist[id]+=1.0;
                if (i==997 && j==999)
                {
                    cout<<part[i].pos[0]<<"  "<<part[i].pos[1]<<"  "<<part[i].pos[2]<<endl;
                    cout<<part[j].pos[0]<<"  "<<part[j].pos[1]<<"  "<<part[j].pos[2]<<endl;
                }
                cout<<i<<"  "<<j<<"  "<<dist<<"  "<<dz<<"  "<<angulo<<"  "<<angulo2<<endl;
                /*Verificamos si se trata de conexion horizontal o vertical
                 * - Horizontal si 150>theta>30
                 * - Vertical en otro caso
                 */
                if(angulo>30.0 && angulo<150.0)
                {
                    *cH+=1.0;
                }
                else
                {
                    *cV+=1.0;
                }
                if(angulo2>30.0 && angulo2<150.0)
                {
                    *cH+=1.0;
                }
                else
                {
                    *cV+=1.0;
                }
                count+=2;
            }
        }
    }
    
    /*Normalizamos el histograma*/
    for(i=0;i<nbins;i++) hist[i]/=double(count);
    /*Normalizamos las conectividades horizontales y verticales con Npart-1
     * Representa el numero de conexiones horizontales y verticales por particula
     */
    *cH/=double(param.npart-1);
    *cV/=double(param.npart-1);
}
#endif	/* CONECTIVIDAD_H */

